Accurate and efficient estimation of local heritability using summary statistics and the linkage disequilibrium matrix

Existing SNP-heritability estimators that leverage summary statistics from genome-wide association studies (GWAS) are much less efficient (i.e., have larger standard errors) than the restricted maximum likelihood (REML) estimators which require access to individual-level data. We introduce a new method for local heritability estimation—Heritability Estimation with high Efficiency using LD and association Summary Statistics (HEELS)—that significantly improves the statistical efficiency of summary-statistics-based heritability estimator and attains comparable statistical efficiency as REML (with a relative statistical efficiency >92%). Moreover, we propose representing the empirical LD matrix as the sum of a low-rank matrix and a banded matrix. We show that this way of modeling the LD can not only reduce the storage and memory cost, but also improve the computational efficiency of heritability estimation. We demonstrate the statistical efficiency of HEELS and the advantages of our proposed LD approximation strategies both in simulations and through empirical analyses of the UK Biobank data.

• On page 9 lines 232-233, the authors mentioned that for fine-tuned approximation settings, the bias in the heritability estimates is even smaller than that when the full LD matrix is used.Are there any explanations why using the approximated LD matrix can have superior performance in terms of bias than using the exact LD matrix?
• The number of SNPs p in Figure 1 is 9,205 and the number of SNPs p in Figure 2 and Table 3 is both 9,220.Since the SNPs in all simulations are from chromosome 22, any explanations for why fewer SNPs were used in the simulations for Figure 1? • All simulations conducted are based on the low-dimensional setting, i.e., the sample size is larger than the number of SNPs.What is the performance of HEELS under the high-dimensional setting?
• On page 15 and the corresponding part in the supplemental notes, the formula for h_SNP^2 is defined as Var(Xi'beta|X)/Var(yi|X).What are Xi and yi? Are they the SNPs and the phenotype for the ith individual, respectively?If so, then the heritability depends on a specific sample, which does not make sense.But based on the formula involving the trace, I am guessing it should be the sum across all individuals in both the numerator and the denominator.Some modifications need to be made here.
• Multiple notation inconsistencies appear in the manuscript, which may cause confusion in reading.
Here are some that I found: 1. On page 16 and the corresponding part in the supplemental notes, when discussing the HEELS procedure, the marginal association statistics and the LD matrix are defined without dividing by root n and n, respectively.These are inconsistent with the ones defined on the previous page.
2. The subscripts and superscripts are used interchangeably in 4 HEELS estimation with unknown sample variance on page 9 of the supplemental notes.The authors mentioned that superscripts are used to distinguish the models for different markers.However, in the formula y=X_j*beta^j+e^j, the subscript was used for the jth column in the SNP matrix X.Such an issue occurred more frequently in the derivation of y'y/n.In the line just above line 159, the sigma_e^j and SE^j should be defined as well.
3. On page 13 in the supplemental notes, how to get from n-beta_J'Rbeta_J to the numerator in the next equation is not clear to me.On page 12, beta_M is defined as R*beta_J.If so, why is there an additional n in the next equation?Maybe it is due to the inconsistency of the notation of the LD matrix mentioned above.Although the authors mentioned beta_J=Sigma^(-1)*beta_M, Sigma is not defined anywhere else.More explanations are needed here.
• There are some typos in the manuscript and in the supplemental notes.Here are some that I found.
1. On page 21 line 514, "…out variants with genotype missingness > 0.01 or have MAF greater than 0.01."Since the main text is talking about common variants, should it be filtering out genotypes having MAF less than 0.01?
Overview: Li et al present HEELS, an approach to estimate SNP heritability using summary statistics and reference LD information.HEELS builds on previous efforts estimating SNP heritability using summary data (e.g., HESS, LDSC, GRE), however rather than use an empirical moment-based solution, it relies on MLE/REML estimates using a score-based approach.Similarly, the manuscript describes a novel approach to represent LD information resulting in large storage savings.The authors perform simulations and real-data analyses to demonstrate the improved statistical efficiency of HEELS over moment-based estimators.I found the manuscript to be exceptionally well written, with helpful balance between details and results split between primary/supplementary info.I have only a few comments.
Major Comments: 1.While the authors are very clear the method assumes only a single genetic variance component (and that extensions are left for future work), it would be helpful to see the performance of HEELS under simulations where MAF, LD, and functional annotations play a role in shaping genetic architecture, and how much information loss there is when inferring a single aggregate parameter.
Response: Thank you very much for your comment.We included additional simulation results to showcase the performance of HEELS compared to GRE and LDSC, when the true genetic architecture is mis-specified.Specifically, for MAF-dependent architecture, we incorporated an  factor to control the strength of MAF dependency and model the SNP-specific genetic variance We observed that consistent with the existing findings of the GREML performance, HEELS is not robust to these mis-specified models and can produce biased heritability estimates (Response Figure 1).These results are not surprising, but instead are aligned with our expectations, given that HEELS behaves almost identically as individual-data based GREML, which has been shown to be biased under mis-specified genetic architectures (Hou et al. 2019, Evans et al. 2018).Note that our simulations relating to the sparse architecture (i.e.twocomponent mixture) can be viewed as incorporating a simple binary annotation in the generative model, but we defer a more detailed analysis of how functional annotations affect our estimator to future studies.Notably, the bias of the HEELS estimator is smaller than that from the LDSC, and the MSE of HEELS is the lowest among the summary-statistics-based methods (Response Figure 2).We have discussed in the Discussion section the future research direction of extending HEELS to allow for the dependence of heritability on annotations.
We have discussed these results in Line 344-349 of the main text and Supplementary Figure 20-21.• First, Benner et al. 2017 focuses on comparing the impact of different LD modeling on fine-mapping whereas our end goal is to ensure the consistency and efficiency of the heritability estimator.Notably, the primary metric evaluated in Benner et al. is the size and the coverage of credible sets of the causal variants.In contrast, we compare the performance of different LD approximation strategies using 1) the accuracy of approximation, defined as the ratio of the norm between the approximated LD matrix and the original empirical LD matrix, and 2) the statistical property of the heritability estimator when the LD is approximated.• Second, we emphasize in our work that the positive semi-definite (PSD) property of the approximated LD has important implications for the HEELS heritability estimator.LDstore does not regard the PSD of the LD approximation as necessary in the shrinkage procedure, and does not consider the factor of PSD in their comparisons.

Response
Please review our changes in Line 185-186 of the main text.
Thank you very much for pointing it out.We have changed the caption.

Response to Reviewer #2:
In this paper, the authors proposed a new estimator, "HEELS", for local heritability.The proposed method attains comparable statistical efficiency as REML by only using summary statistics.Moreover, the authors also proposed a unified framework to approximate the empirical LD matrix by decomposing it into the sum of a banded matrix and a low-rank matrix.Such an approximation can not only reduce the storage and memory cost of using the LD matrix, but also improves the computation efficiency of the proposed HEELS.Overall, the methodology is sound and the paper is well-written.
Response: Thank you for your positive comments on our paper.
Below are my comments and concerns on the paper.
• On page 9 lines 232-233, the authors mentioned that for fine-tuned approximation settings, the bias in the heritability estimates is even smaller than that when the full LD matrix is used.Are there any explanations why using the approximated LD matrix can have superior performance in terms of bias than using the exact LD matrix?
Response: Thank you very much for your comment.The improved accuracy of the heritability estimator under certain settings when LD is approximated may be attributable to regularization, where the approximation leads to an estimate of the LD that is closer to the true covariance.
Intuitively, the regularized estimators will perform better if the true covariance matrix has structures that are well-captured by constraints or assumptions imposed by the regularization, such as sparsity or banded-ness.Seminal works on covariance matrix estimation have provided high-probability analyses that bound the errors of regularized estimators when the true covariance is sparse or banded (Section 6.5 of Wainwright 2019, Cai et al. 2010, Bien et al. 2016).
In our case, we hypothesize that the LD approximation strategies that lead to more accurate heritability estimates might have better matched the original LD (in particular, the specific banded + low-rank structure assumed).In a recently published paper (Salehi et al. 2022), researchers have found that LD regularization can lead to less biased BLUP estimates.Although this paper examines the impact of the regularized LD precision matrix on the SNPeffect predictors, the better performance of the regularized LD precision matrix is consistent with our observation of lower bias in the heritability estimates in the presence of regularization.
• The number of SNPs p in Figure 1 is 9,205 and the number of SNPs p in Figure 2 and Table 3 is both 9,220.Since the SNPs in all simulations are from chromosome 22, any explanations for why fewer SNPs were used in the simulations for Figure 1?
Response: Thank you for pointing it out.The number of SNPs is different between these two sets of analyses because we applied the MAF > 0.01 filter separately to two different samples: • The results in Figure 1 are generated based on a random subsample of the full UKB sample (n=30,000).We restricted analyses to the smaller sample due to the limited scalability of GCTA when applied to biobank sized datasets.• On the other hand, the results in Figure 2 and Table 3 are based on using the full UKB sample (n=332,430).

Figure 1 .
Distribution of the heritability estimates from HEELS under mis-specified models.The simulation results are based on real genotypic data from random subsets of unrelated individuals in the UK Biobank, array SNPs on chromosome 22 with MAF > 0.01.The  value determines the strength of the MAF-dependency, with   2 ∝ (  (1 −   )) 1+ .Setting  to -1 corresponds to the GCTA model; setting  to -0.25 corresponds to the LDAK model.LDdependent architecture: assume the genetic variance of a variant is inversely proportional to its level of LD.Red-dotted line: true heritability of 0.25.The row panels represent polygenicity or the proportion of SNPs with non-zero effects.Response Figure 2. Comparison of MSE between different summary-statistics-based heritability estimation methods (HEELS, GRE and LDSC) under mis-specified models.The simulation results are based on real genotypic data from random subsets of unrelated individuals in the UK Biobank, array SNPs on chromosome 22 with MAF > 0.01.The  value determines the strength of the MAF-dependency, with   2 ∝ (  (1 −   )) 1+ .Setting  to -1 corresponds to the GCTA model; setting  to -0.25 corresponds to the LDAK model.LDdependent architecture: assume the genetic variance of a variant is inversely proportional to its level of LD.Red-dotted line: true heritability of 0.25.The row panels represent polygenicity or the proportion of SNPs with non-zero effects, e.g., 0.1 means 10% of the SNPs are causal.